# a function that performs the Lilliefors test for normality
# version: based on Abdi & Molin (2007) approximation of the Lilliefors distribution
# source: H. Abdi & P. Molin: Lilliefors / Van Soest Normality Test - https://personal.utdallas.edu/~herve/Abdi-Lillie2007-pretty.pdf
# see also: "Lilliefors Test with Abdi-Molin Approximation" by Masayuki Tanaka - https://www.mathworks.com/matlabcentral/fileexchange/57832-lilliefors-test-with-abdi-molin-approximation
# ENSURE that Octave can find the function, use the "addpath" command: addpath ("/[insert the path to your Octave folder]/Octave")
# NOTE: to run this function you may also need to load the statistics package, "pkg load statistics"

# function arguments:
# x - a vector of data for which the normality assumption is tested
# function outputs: pval - p-value, Dn - Kolmogorov-Smirnov statistics, unscaled
# NOTE: the interpretation is as follows: if p-value > alpha, the null hypothesis - that x is normally distributed - cannot be rejected, otherwise it should be rejected

function [pval, Dn] = Ltest2 (x)
if (isvector(x) == false)
    error ("the input must be a vector")
endif
N = length(x);
if (isrow(x) == true)
    i = 1:N;
elseif (iscolumn(x) == true)
    i = transpose(1:N);
endif
s = sort(x);
fx = normcdf(zscore(s));
D1 = max(abs(fx-i/N));
D0 = max(abs(fx-(i-1)/N));
Dn = max(D1,D0);

# The first step: compute a quantity  A
b0 = 0.37872256037043;
b1 = 1.30748185078790;
b2 = 0.08861783849346;
A = (-(b1 + N) + sqrt((b1 + N)^2 - 4*b2*(b0 - Dn^(-2)))) / (2*b2);

# The second step: compute a polynomial that estimates the probability associated to a given value of Dn
c = [0.00000575586834, -0.00011872227037, 0.00069650013110, 0.00287462087623, -0.06353440854207, 0.40466213484419, -1.39874347510845, 2.80015798142101, -3.02959249450445, 1.67819837908004, -0.37782822932809];
pval = polyval(c, A);

if (pval < 0)
    pval = 0;
endif
if(pval > 1)
    pval = 1;
endif

endfunction
